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the lines of the well-known inf criterion for the factorization method. 
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1 Introduction 




Figure 1: Reference setting (left) and actual setting (right). 

We consider an inverse scattering problem consisting in shape reconstruction from physical measure- 
ments. These problems are generally non-linear and ill-posed. More specifically we address the problem of 
reconstructing the support of a perturbation in a given inhomogeneous background medium from acoustic 
far-field measurements generated with plane waves. It may indeed happen that, in some places, the actual 
index is different from the reference value, as seen in Figure 1. This could happen for instance from a 
deterioration or an incorrect estimation of the actual index. We then say that there is a defect at any 
point where the reference index is different from the actual index. 

A wide range of methods achieve the localization of obstacles by a sampling approach: the points of 
the unknown domain are characterized by a binary test that has to be applied to the whole space. For 
most of them, the first formulation of this pointwise test is to check if some well chosen test-function is 
in the range of an implicitly defined operator. See [4, 15] and references therein for a topical review. A 
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natural way to proceed is then to connect the range of the implicit operator to the range of an operator 
explicitly defined from the actually available measurements. This is the principle of the linear sampling 
method [5, 7, 3] or of the factorization method [9, 12, 1]. Yet, when looking for perturbations in non- 
homogeneous background media, it is only recently that a factorization method has been proposed to 
reconstruct the shape of defects [14, 8]. 

However, we investigate in this paper an optimization approach in the lines of the inf criterion [11] to 
deal with the implicit operator's range. We show that this leads to a characterization of the defects as the 
support of the following function: 

M w (z) := inf (/„,(*), <J> G L 2 (S d ~ r ) and (*, u no (-,z)) = l) , 

where u no is the total field for the reference index and the form fw is given by 

f w (9) := \(W% *) ia(si - 1) 

with a measurements operator W explicitly built from the available data. Throughout this paper, (•, -) x 
stands for the usual hermitian inner product in the Hilbert space H. 

This paper is structured as follows. In section 2, the mathematical setting is specified and the implicit 
localization of the defects is recalled. This localization is then expressed as a binary pointwise test involving 
an constrained optimization problem in section 3. Finally, in section 4, we investigate numerical methods 
for solving this optimization problem. We end by some conclusions. 



2 Formal localization 

If we consider time-harmonic acoustic waves with a fixed wave number fc, the spatial part of the wave 
equation is modeled by the Helmholtz equation [6]. Inhomogeneous media are then represented by an 
acoustic refraction index, denoted by n G L°°(R d ), and so the total field, denoted by u n , is assumed to 
satisfy 

Au n + k 2 n(x)u n =0, x G R d , (1) 

where d is the problem's dimension (d = 2 or 3). We consider compactly supported inhomogeneities and 
denote by D the support of n(x) — 1. Also, we denote an incoming wave satisfying (1) with n = 1 by 
u % G L 2 oc (K d ). The total field is then the sum of this incoming wave and the wave scattered by the 
inhomogeneous medium, denoted by u s G L 2 oc (R d ): 

u n := u s +u\ (2) 

where the scattered wave is assumed to satisfy the Sommerfeld radiation condition 

d r u s =iku s + o(\x\~^^j . (3) 

Then, the linear system (l)-(3) defines u n uniquely from u l , and it is known to be invertible in L 2 (D). 
Thus, let us denote the corresponding automorphism by 

T n : L\D) -> L 2 (D), 

Besides, the outgoing part of a wave has an asymptotic behaviour called the far field pattern, denoted 
by G L 2 (6' (i_1 ), see figure 2, and given by the Atkinson expansion [17] 

u n (x) —u^x) +7-^^<°(£) + o(Vr^ 1 ) , x:= ^ eS d ~\ 

\x\— V 1 N 
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Figure 2: General setting and notations. 



where 7 depends only on the dimension and is defined by 

7 := 



if d = 2, 



8-Kk 

h if d = 3. 
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Furthermore, for practical reasons, we will mainly consider scattered waves having a plane-wave source. 
These plane- waves are defined by 

u l (0,x) := exp(ikO • x), 

where is a unitary vector in S^ 1 as depicted on figure 2. Then, let us denote the total field at the point 
x G M d and with a plane- wave source of incoming direction by 

Un(0,x) :=Tn(u\9,-))(x). 

The corresponding far-field pattern in the measurement direction x G S^ 1 will be denoted by u™(0,x). 
Lastly, these measurements will be used in the form of the classical far-field operator F n : L 2 (6' d_1 ) — >• 
L 2 ^- 1 ), defined by 

F n g(x) := (g, u™(-,x)) . 

We want to reconstruct the shape of defects in a reference medium whose index is denoted by no G 
L°°(D). Let then ni G L°°(D) denote the actual index, altered by the presence of these defects. So, 
denote the support of the difference between the two indices (see figure 1) by 

Q := support (774 — no). 

Thus, the goal is to reconstruct the domain through its associated characteristic function denoted 1^, 
from the reference index no and far-field measurements . Yet, it has been shown that the unknown 
domain can be characterized by a set of test functions and the range of some implicitly defined operator. 

Theorem 1 (Implicit domain characterization). [8, Theorem 3.2] Let us define the operator C : L 2 (D) — » 
L 2 ^- 1 ) by 

Cf(x) = (f,U no (x,-)) L 2( D y 

Then, for each z G we have 

zen un (.,z) eft(Cln), (4) 
where 1q : L 2 (Q) — )► L 2 (D) is the restriction to Q given by 



0, x£Q 
f{x), x G VL 
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3 Explicit identification of the defects 



This section gives an explicit formulation of characterization (4). To do so, we proceed in three steps. First, 
we introduce a practical result formulating the belonging to the range of any given (bounded) operator L 
as an inf criterion based on any function related to L through specific assumptions. Then, we construct 
such a well-suited function for the operator we are interested in, namely CIq. Finally, we can give the 
explicit characterization of the defect Q. 



3.1 Formulation of a bounded operator's range through an inf criterion 

To deal with the range of Cl^, depending on the unknown domain ^, let us recall the following charac- 
terization of an operator's range. 

Lemma 2 (Range characterization). [14, Lemma 2.1] Let L : Hi — >> H 2 be a bounded operator between 
two Hilbert spaces H\ and H 2 , and let cp G H2. Then, noting H 2 the dual of H 2 and identifying Hi with 
its dual, (p G 7£ (L) if and only if there exists c > such that for all \£ G H 2 

|<*, ^>|<c||L**||. (5) 

Hence, any form comparable to ||L*(-)|| can be used to characterize the range of the operator L. 

Corollary 3. With L as in lemma 2, let f : H 2 ^ be comparable to L* in the sense that there exists 
ci > and c 2 > such that 

ci ||L*tf || < /(*) < c 2 ||L**|| , * e HI (6) 

Then, ip G 1Z (L) if and only if there exists C3 > such that for all ^ G H 2 

|<tf, <p)\ < c 3 /(tf). (7) 

Proof The result (7) is a straightforward combination of the characterization (5) and (6). □ Finally, 
from corollary 3, we deduce the following characterization based on an inf criterion. 

Corollary 4. With f as in corollary 3, then cp G 7Z (L) if and only if 

< inf {/(*), * G H$ and (*, ip) = 1} . 

Proof First, if (p 7£(L), from corollary 3, for each c > there exists Sk c G H 2 such that |(\P C > tp)\ > 
c/(* c ). Since f(SSf c ) ^ by (6), then (^ c , <p) ^ and we can define ip c := ^ c /(^ c , ip). From (6), it follows 

So, there is a set of functions ijj c G H 2 , satisfying (^ c , y?) = 1, such that f(i/j c ) ~~ ^ when c — » 00. 

Next, let (p G H 2 \ {0} be in the range of the operator L and let ^ G #2 satisfy = 1. Thus, from 

corollary 3, there exists C3 > such that 1 ^ csf(^) and the infimum is not vanishing. Finally, if cp = 0, 
which is always in the range of L, the infimum is evaluated over an empty set and will then conventionally 
be given the value +00. □ 
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3.2 Construction of a well-suited objective function 



As shown in corollary 4, the range of a linear bounded operator L can be formulated as a usual constrained 
optimization problem without the exact knowledge of L. Indeed, it suffices to find any form / satisfying (6). 
Hence, to get an explicit characterization of the domain Q from theorem 1, we look for a form / satisfying 

ci||l n C*tt|| </(*) < ca ||1„C**|| , ta 2 ^" 1 ). (8) 

To achieve this, following [11], we consider the form fw ' L 2 ^^ 1 ) —> M + defined by 



f w (9) := (W9, * } i2 



(9) 



where W is a measurements operator explicitly built from the available data. 

A first guess for the operator W would be to consider (F ni — F no ): the difference between the far- field 
operators corresponding respectively to the reference index no and the actual index n\. But it has been 
shown in [8] that this subtraction does not yield some crucial factorization. Thus, we have to restrain 
ourselves to full bi-static data (u™ and are known over S^ -1 x S^ -1 ) and non- absorbing media (tiq(x) 
and ni(x) G M) so we can consider the operator W : L 2 (5' d_1 ) — » L 2 (S' d_1 ) defined by 



W 



:= (id + 2^| 7 | 2 F no )(F ni -F no ). 



Under these assumptions, it has been shown that this measurements operator W has a factorization of 
the form (see [8, Corollary 4.7, Corollary 4.3]) 



W = CAC\ 

where the operator A is an automorphism on L 2 (Q) defined by 

A := l^/c 2 (ni - riQ)T ni T~ 1 ln. 

So, we have 



(10) 



(11) 



(12) 



and as such, the inequalities (8) are then the continuity and the coercivity of the operator A, at least on 
the range of C* . We are now going to show that this coercivity is related to the contrast between the 
reference index no and the actual values of ni, i.e. the defects should be clearly distinguished from the 
background. We thus make the following geometrical assumption. 

Assumption 5. Assume that no and rt\ are real valued and that either (ni — no) or (no — ni) is locally 
bounded from below: 

• for any compact subset uo included in there exists c > such that (711(2) — rio(z)) ^ c for almost 
all z G uj, 



• for any compact subset uj included in f], there exists c > such that (no (2) — ni(z)) ^ c for almost 
all z G uj. 

Moreover, for a fixed geometry of defects, some wave numbers k may produce resonances that cancel the 
outgoing wave. Indeed, the operator mapping an incoming wave to the corresponding far-field pattern is 
not one-to-one in the case of inhomogeneous media. This corresponds to the so-called interior transmission 
eigenvalues arise [10]. 
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Definition 6. We call k an interior transmission eigenvalue for the pair of indices (no, n\) if there exists 
a non-vanishing (source, solution) pair, denoted by (/i, u) G (L 2 (£l)) , such that 

' (A + k 2 no)u — —k 2 (ni — no)h in 
(A + k 2 n!)h = in ft, 
m = on <9ft, 
d v u = on dfl. 

Since we want to avoid these "pathological" values, it is useful to know that this is a rare case. 

Lemma 7. // the indices no and n± are real-valued, then the set of interior transmission eigenvalue for 
the pair of indices (no,ni) is discrete. Furthermore, if there are infinitely many, they only accumulate at 
+oo. 

Proof The proof follows exactly the lines of [13, Theorems 4.13 and 4.14], by adapting the notations. □ 
With these two geometrical restrictions, the coercivity of the operator A on the range of C* is then obtained 
in lemma 9, using the following result. 

Lemma 8. [13, Lemma 1.17] Let Y be a subset of a reflexive Banach space X and A, Aq : X* — >• X be 
linear and bounded operators such that 

1. ((/?, Acp) G C \ (—oo, 0] for all cp ^ in the closure ofY 

2. ((/?, A^cp) G R and there exists Co > with ((/?, A cp) ^ cq \\(p\\ for all (pinY 

3. A — Aq is compact. 

Then, there exists c > such that for all cp G F it holds \(A(p, cp)\ ^ c • 

Finally, we group the properties of operator A in the following lemma. 

Lemma 9. Under assumption 5, ifk is not an interior transmission eigenvalue in the sense of definition 6, 
then the operator A, defined by (11), is coercive on the range of C* , defined in theorem 1. Namely, there 
exists c\ > such that 



ci||l n C7**|| La(n) ^ 



(AC*V, C**) L 2 m i G L (S ). (13) 



Proof. This is a straightforward application of lemma 8 with Y = IZiC*) and X — L 2 (Q). The required 
assumptions have been partially shown [8, lemma 5.3] but, for convenience, we give a complete proof. 

1. Choose cp G H(C*). By definition of operator C, this is a total field for the refraction index no- 
Hence, there exists an incident field / such that cp = T no (f). Let us set u no = cp and u ni = T ni (f). 
Thus, we obtain Acp = k 2 (ni — no)u ni . Moreover, choosing R > such that the ball Br of radius R 
contains it holds that 

/ & 2 (ni - n )u ni (u no - u ni ) 
Jn 



(A + k n )(u ni - u no )(u ni - u no ) 



I k 2 n \u ni - u no \ 2 - \V(u ni - u no )\ 2 + / (u ni - u no )d u (u ni - u no ). 

J Br J Sr 
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By letting R go to infinity, it comes 

k 2 (ni - n )u ni (u no - u ni ) 



k 2 n \u ni -u no \ 2 -\V(u ni -u no )\ 2 +ik\j\ 2 I \u™ -<°' 2 



m ^n 1 

S d-1 



Hence, taking the imaginary part yields 



This shows that Im(A(p,(p) ^ and if this quantity is vanishing, we deduce that = . 
Moreover, out of we have 

(A + k 2 n )u no = (A + /c 2 ni)ix ni = (A + k 2 n )u ni . 

As a consequence, the unique continuation principle [6, theorem 8.6] yields u ni = u no out of Q. So, the 
quantity w := (u ni — u no ) has its support included in Q and satisfies (A + k 2 no)w = —{n\ — n$)u ni . 
If k is not a transmission eigenvalue, we then have w = and u ni |^ =0. Finally, this implies 
that u ni = by the unique continuation principle and all these results are extended to 7£(C*) by 
continuity to prove item 1. 

2. Furthermore, we also have that 

(Aip,<p)= / k 2 (ni - n )|^ no | 2 + / k 2 (nr - n )(u ni - u no )u^ 
Jn Jn 

with Aq = k 2 (ni — no) I and K = k 2 (ni — no)(T ni T~ 1 — I). Under assumption 5, Aq is clearly 
coercive and self- adjoint. 

3. Moreover, K = k 2 (ni — no)(T ni — T no )T~^ . Yet, it is known from the Lippmann-Schwinger equa- 
tion [6, equation (8.12)] that T n = id — TT ni where T is some compact operator. Thus, (T ni — T no ) 
is compact, and so is K. 

□ 



3.3 Characterization of the defects from the measurements 

We are now able to state an explicit localization of the defects in the form of a constrained optimization 
problem that extends the characterization proposed in [11]. 

Theorem 10. Assume that k is not an interior transmission eigenvalue for the indices no et n\, following 
definition 6, that these indices are contrasted following assumption 5 and that we have full bi-static data 
(i.e. u^ Q and are known over S^ 1 x S^ 1 ). 
We can then define the value 

M w (z) :=inf(jW03>), * S L 2 (S d ~ r ) and (*, u no (; zj) = l) , (14) 

[ \ /L 2 (S' d - 1 ) J 

and for each point z G M d we have 

z e M w (z) > 0. (15) 
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Proof. Theorem 1 characterizes the domain of the defects Q by 



zeQ ^ u no (-,z) £ 1Z (CIq) . 

The result (15) is then a direct consequence of corollary 4 with / = fw : <p = u no (-,z), L = CIq and 
B.2 = L 2 (6' (i_1 ). So what is left to prove is that the double inequality (8) holds. Now, we deduce from (12) 
that the right inequality comes from the boundedness of the operator A. The left one was established in 
lemma 9. □ 



4 Numerical methods for the computation of the infimum's value 

In theorem 10, we have expressed the localization of the defects as a pointwise binary test, taking the 
form of a constrained optimization problem. We are now interested in the numerical computation of the 
values of the function Mw( z ) defined in (14). We thus explicit, and then test, two usual minimization 
algorithms working on this problem: the steepest descent and the gradient projection. 



4.1 Algorithms 

Both minimization methods we are going to use require explicit gradients of the cost function. To simplify 
their expression, we will consider the minimization of the following form on L 2 (5' d_1 ), defined for any \£ 

by 

:= (jV(*)) 4 = \(WV, *> i2{s «.-i) 12 

Since we only want to know if the infimum is vanishing, this gives results equivalent to (14) and we thus 
have to evaluate 



(M w (z)) 4 := inf {/^(*), * e L 2 (S d ~ 1 ) and (*, ^o(^)) L2(sd _ 1} = l} 



(16) 



Remark 11. The function Mw{z) has been proved to vanish outside of the defects and yet, it can be 
seen that the value can never be attained while satisfying the constraint. 

To show this, let £ L 2 (6' (i_1 ) be such that f^i®*) = 0. Then, from factorization (10) and inequali- 
ties (13), it holds that /^(#*) = (AC*®*, C*^*) 2 ^ ||1^C*^^|| 4 , that is \n = 0. Furthermore, it 
is easy to see that C*\I/* satisfies the Helmholtz equation (1) on M d with n = n$. The unique continuation 
principle then yields C*\I/* — 0, and thus = by inject ivity of C*, which has been shown in the proof 



of [8, Proposition 5.4]. As a consequence, can not satisfy u no (-,z)^ = 1. 

This points out that numerical approximations of a vanishing inf might be mistaken with exact non-zero 
inf values. Some care will thus have to be taken to set them apart, so that the plot of M\y(z) can be 
used to localize the defects. 

We now turn to the feasible set. For z £ R d , let us denote by C z the affine hyperplane 



1 hn (;z)\\ X /LHS*-i) J 

The orthogonal projection on C z , denoted by Pc z1 is defined on L 2 (6' d_1 ) by the affine mapping 

Pc z ^ := * - ^ /% (-, z) + u ^> z > (17) 

hn (-iZ)\\ IKo('^)ll 
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Hence, finding an infimum of over C z is equivalent to looking for the infimum over L 2 (S d x ) of 

We can therefore compute a minimizing sequence x n for the form P^ by any unconstrained optimization 
method. For convenience, all gradients and hessians are calculated in section A.l and their finite dimension 
formulation is given in section A. 2. It then follows from theorem 10 that if {Pf^( x n)) n goes to 0, the 
point z is outside ^, and inside otherwise. A basic example of descent is presented in algorithm 1: 

Input: xq chosen in L 2 (5' (i_1 ) 
repeat 

Compute a n such that P^(x n — a n V Pf^(x n )) < P^{x n )] 

Update x n+ i <- x n - a n VP f ^{x n )] 
until ||x n+ i - x n \\ /(l + ||x n ||) < e\ 
Output: P f 4 v (x N ) 

Algorithm 1: Steepest descent 
Moreover, since the projection on C z is easy to write, we can also consider a gradient projection 
method [16]. As previously, if (fw( x n)) n g° es to 0, the point z is outside Q. The principle of the 
gradient projection is presented in algorithm 2: 
Input: xo chosen in C z 
repeat 

Compute a n such that f^(x n - a n Vf^(x n )) < f^(x n ); 

Project and update x n+ i <- Pc z ( x n ~ a n \/f^(x n )); 
until \\x n +i - x n \\/(l + ||x n ||) < e\ 
Output: fw(xN) 

Algorithm 2: Gradient projection 

Remark 12. Since the projector Pq z is an affine map, the proposed steepest descent and gradient pro- 
jection methods are very close. Indeed, we note that by choosing a constant descent step a n = a and a 
starting point xq g C z , both algorithms define the same sequence. The computation of a n is however done 
before the projection in algorithm 2, and after it in algorithm 1. Thus, the values coming up for a n in 
each of the proposed algorithms will seemingly be different and produce different sequences x n . 

4.2 Numerical validation on a simple case 




Figure 3: A simple study case 

We present here some 2D numerical results in the simple case illustrated on figure 3. The considered 
object of support is a disc D of section 2.1 containing defects which are also in the shape of a disc £1 of 
section 0.6. With a fixed wave number k = 10, the size of the object is then thrice the wavelength and 
the size of the defects is approximatively one wavelength. Also, the reference index no takes its values in 
[1.56, 1.84] inside D and the perturbed version m takes its values in [2.01, 2.16] inside ft. Finally, we used 
99 incoming/measurement directions evenly distributed over [0, 2tt]. 
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(a) Steepest descent algorithm (b) Gradient projection algorithm 




(c) Matlab's fminunc function (d) Matlab's fmincon function 



Number of iterations 


min med max 400) 


Steepest descent 
Gradient projection 
Matlab's fminunc function 
Matlab's fmincon function 


11 400 400 
14 120 400 

8 17 68 

9 19 58 



Figure 4: Values of log 10 M\y(z) computed by various optimization methods with a relative tolerance on 
x n set to 10 -9 as the only stopping rule 
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Figure 4 displays the infimums's values [M.w{zj) in a log 10 scale, respectively obtained for each 
sampling point Z{ (about 8500) with algorithm 1 (figure 4a) and algorithm 2 (figure 4b). Note that these 
sampling points are unrelated to the finite elements nodes that were used to generate the data. The 
linear search for the step length a was done by simple dichotomy on the gradient. We also present the 
results obtained with Matlab's f minunc function applied to the form (figure 4c) and Matlab's f mincon 
function applied to the form over the feasible set C Zi (figure 4d). The f minunc and f mincon functions 
are based on the interior-reflective Newton method described in [2]. Furthermore, for each sampling point 
and each algorithm, the sequence has been initialized by xq = Pc z .(0) = u no (-, Zi)/\\u no (', Zi)\\ 2 G C Zi . 
Indeed, it seems natural to start with a point already satisfying the constraint. 

As can be expected, Matlab's functions give much faster results, but we note that even the very basic 
algorithms we proposed yield acceptable results. This shows that the iterative optimization approach 
resulting in theorem 10 can produce a satisfactory localization of the defects, but for a computational cost 
that is not controlled at this point. 

To go further, we also tested the sensitivity with respect to the data, to take into account simulation 
or measurement inaccuracy. This was done by adding uniform random noise to the measurements and 
thus, using u e n such that \u e n — || < e || izj£ ||. The optimization approach turns out to be quite robust 
regarding this criterion. Indeed, figure 5 shows that acceptable results are still obtained with 10% relative 
noise, with a seemingly better visualization for the gradient projection algorithm. As this is sometimes the 
case, we also note that adding some noise has a slightly regularizing effect that visibly enhances convergence 
speed for the basic algorithms 1 and 2. 

Remark 13. The results in figures 4 and 5 were obtained by considering only the relative variation on 
x n , as proposed in algorithm 1: \\x n +i — x n \\ /(l + ||x n ||) < 10 -9 . Yet, many optimization algorithms rely 
on multiple stopping criterions, including a relative variation tolerance with respect to the cost function's 
values. We see here that these values are very small and thus hardly usable in a stopping rule. The same 
goes for the gradient and first order optimality criteria. As a consequence, we had to set the relative 
tolerance regarding the cost function to 10 -32 in Matlab's optimization functions, so that this condition 
would never be triggered. With the tolerances for the relative variation regarding x n and the cost function 
set to their default (10 -6 , see Matlab's help), we see in figure 6 that Matlab's functions produce the 
opposite result to what was expected. Indeed, we see that the values of Mw{z) outside ft are close to zero 
but higher than the ones inside Q. On the other hand, the simple algorithms presented in the previous 
section still yield the expected localizations, and at a lower computational time. This highlights that a 
special care has to be taken regarding the stopping rule, as we compare the minimal values issued from 
different minimization problems. 

4.3 Experimentations on a non-trivial absorbing example 

It is illustrated in figure 7e that comparable results are obtained on a more elaborate example with two 
non-convex and non-connected defects depicted in figures 7a-7d. 

Besides, theorem 10 is stated under some physical restrictions arising from the use of the scattering 
operator (id + 2ik \j\ 2 F no ). The numerical methods proposed in this section can however straightforwardly 
be extended to absorbing media and limited far-field data. We see in figure 7e that defects in complex 
valued indices are still correctly localized, even with 2% uniform random noise added to the measurements. 

Finally, even with limited far-field data, the gradient projection algorithm provided a satisfactory recon- 
struction, displayed in figure 8a. In this last example, the 99 evenly distributed incidence/measurement 
directions were taken in [0, |7r]. Still, even with this fairly high amount of data, it is to be noted that 
the four algorithms presented in this paper do not yield comparable results in this case. Indeed, we see 
in figure 8b that Matlab's functions fail to provide a usable reconstruction, despite our testing on a wide 
range of optimization parameters. 
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(a) Steepest descent algorithm (b) Gradient projection algorithm 




(c) Matlab's fminunc function (d) Matlab's fmincon function 



Number of iterations 


min med max 400) 


Steepest descent 
Gradient projection 
Matlab's fminunc function 
Matlab's fmincon function 


5 10 400 

6 12 145 
6 10 14 
8 12 18 



Figure 5: Values of log 10 Aiw( z ) with 10% uniform noise added to the measurements 
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(a) Steepest descent algorithm 



-2-10 1 

(b) Gradient projection algorithm 




-10 12 -2-10 1 



(c) Matlab's fminunc function (d) Matlab's fmincon function 



Number of iterations 


min med max 400) 


Steepest descent 
Gradient projection 
Matlab's fminunc function 
Matlab's fmincon function 


4 46 400 

5 29 400 
4 9 23 
7 9 26 



Figure 6: Localization of the defects with the relative tolerances set to 10 
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(a) Reference index's real part 
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(b) Reference index's imaginary part 
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(c) Perturbed index's real part 
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Figure 7: Reconstruction of a non-trivial perturbation with the gradient projection algorithm 
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(a) Gradient projection algorithm (b) Matlab's fminunc function 
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Figure 8: Reconstruction of a non-trivial perturbation with 2% and limited far-field data 

5 Conclusion 

We have characterized the localization of defects in an inhomogeneous reference index by an optimization 
problem that is built only on the available data. Objective function and feasible set turn out to be very 
simple. This problem can thus be solved through a wide range of well known optimization methods, 
which we have numerically illustrated four examples of. Yet, some limitations were noticed, regarding the 
convergence speed on some cases and the required amount of data. Issues for which successful results were 
obtained using spectral methods in [8]. This opens the perspective of looking for some stabilization, or 
more robust versions of the proposed optimization algorithms. 



Appendix 

A.l Derivatives of the objective function 

We give here the derivatives of the functions involved in the computation of Mw(z) as defined by (16). However, 
since the values of the form fw are real, the differential can not be C-linear. We therefore have to split the elements 
of L 2 (S d ~ 1 ) into their real and imaginary parts and consider the objective function on pairs of real- valued functions 
to obtain proper R-linear differentials. The induced gradient is then given in the following lemma. 

Lemma 14. The gradient of the form : L 2 (5 d_1 ,R) x L 2 (S d ~\R) -> R, defined by 
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is given by 



:=4 



V 



Re Gf 
Tm Gf 



(A-l) 



where the function Gf is an endomorphism on (L 2 (S d 1 ,M)) 2 defined, with help of the self-adjoint parts Wr := 
(W + W*)/2 andWi :=(W -W*)/2i, by 



(A-2) 



A/so ; £/ie second derivative of f w defined on (L (S ,R)) , /or eac/i pom* (</>,^) € (L (5 , R)) , 6j/ 



= 4 



Re 
Im 



a\\ 



b 

a 

bJJ 



(A-3) 



where the operator 



is an endomorphism on L (S ) defined for each [ ^) e (L (S )) 6?/ 



2Re (VKr^, $) WrV + (W^, ) + 2Re (WiV, $) + (W/tf , ) Wj$, 



where \I> := (f> + and $ := a + i&. 



Proof. First, for (J J and (J J in L 2 (5 d_1 ,R) x L 2 (5 d_1 ,R) we denote 
Let then : L 2 (5 d_1 ,R) x L 2 (5 d_1 ,R) ^ C be denned by 



so we have \jpj ~ ^ w \^<0y ' Moreover, the operator is not self-adjoint but is nevertheless an endomor- 
phism. Hence, we can use its real and imaginary parts as defined by 

Wr := (W + W*)/2, Wi := (W - W*)/2i. 



It follows that (WrV, V) = Re (WV, tf) and (WiV, tf) = Im Since the operators et Wi are 

self-adjoint, with & := u + iv we obtain 



<t>\ (u 



2Re (W^, $) + 2iRe {WiV, $) . 



Hence, the differential of is given by 



<A (u 



= 2Re ((WrV, ) + i (Wj*, ) , 2Re (Wj^, $) + 2iRe (Wj*, $)) 

- 4Re ({WrV, tf) W^tf + (W/tf , tf) W/*, $) 
= 4Re (G/(tf), 
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where Gf is defined by (A-2). The gradient (A-l) is then written by recalling that for hermitian inner products, it 
holds that 

-</.H(wMS)>' 

Finally, we get the second derivative by differentiating the gradient. Since the operators Wr and Wi are C- linear 
and self-adjoint, it comes 

DG f (V)($) = 2Re (WrV, $) WrV + (WrV, W R ^ + 2Re (WiV, $) WiV + (W/tf , 

□ 

As a consequence, we also have to adapt the projection Pc z (17) and use its counterpart Pc z defined as an affine 
endomorphism of L 2 (S d ~\R) x L 2 (S d ~\R) by 

-* + * (A - 5) 

The gradient of the objective function's projection, denoted by Pf4^ , is then given by 

v Pf4 , TI = ,^(p c .rG/(ft.*) v = 

' Im (P C , )*G/(P C , *] ' 



where Pc z is the linear part of the projection and where Gf is given by (A-2). 

A. 2 Finite dimension approximation 

For the finite dimension approximation, the operators Wr and Wi have a complex matrix representation. In 
order to write the gradient in more natural terms of matrix- vector products, we thus denote Wr and Wi the 
corresponding real valued expanded matrices defined by 

Wr _ /Re (W R ) -Im (W R )\ ^ ^ _ f Re (Wj) -Im (Wi) 



Im (Wr) Re (Wr) J ' 1 " V Im ( W i) Re (W> 

Moreover, we assume the standard change of basis M2, where M is a discretization of the inner product, so that 
the transposition correctly represents the adjoint of operators. Furthermore, with and ip two elements of the 
discretized version of L 2 (5' <i_1 ), in the basis M^ 5 denote 



With these notations it comes that 



Wrx = 'Re W R ($ + #)\ ^ Wix _ fRe + 



Im W R ((j) + iijj)) ' 1 VlmW>(0 + #) 



and by recalling (A-4), we have 



x T W R x = Re (Wr{(/> + #), (0 + #)) = {Wr{4> + #), (0 + #)> , 
x T Wix - Re (Wj(0 + #), (0 + #)) = (Wj(0 + #), (0 + #)) . 

It follows from (A-l) that the (expanded) gradient of the form is given by 

Vfw(x) = 4W R x(x T W R x) + 4Wix(x t Wix). 

It also follows from (A-3) that the corresponding hessian matrix (x) is given by 

H /4 (x) = 8W R x <g> Wrx + 4W R (x T W R x) + 8W T x <g> Wix + 4Wi(x t Wix) 
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/6i 




and b := 




{'J 







where the tensor product between two column vectors a := a2 and b := 2 is defined by the matrix given 

W w 

in columns by 

a®b:= (a6i afo . . .) . 

Finally, it then comes from a straightforward calculation that the finite dimension approximations of the gradient 
and the hessian matrix for the form P f 4 are 

VP f 4 ? (x) = p T Vf^(Pc,x), 



and 



— >-T 



Hp (x) = P H^(P Cz x)P, 

J w J w 

where P is the expanded matrix representation of the projection's linear part. 
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